Number-conserving master equation theory for a dilute Bose-Einstein condensate 
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We describe the transition of N weakly interacting atoms into a Bose-Einstein condensate within 
a number-conserving quantum master equation theory. Based on the separation of time scales 
for condensate formation and non-condensate thermalization, we derive a master equation for the 
condensate subsystem in the presence of the non-condensate environment under the inclusion of all 
two body interaction processes. We numerically monitor the condensate particle number distribution 
during condensate formation, and derive a condition under which the unique equilibrium steady state 
of a dilute, weakly interacting Bose-Einstein condensate is given by a Gibbs-Boltzmann thermal state 
of N non-interacting atoms. 
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I. INTRODUCTION 

After almost one century of theoretical works to under- 
stand the existence, analysis and creation of a new state 
of matter at ultracold temperatures, the first experimen- 
tal observation of a Bose-Einstein condensate was pre- 
sented in Refs. [3, Hj]. Nowadays, Bose-Einstein conden- 
sates are well established as a distinguished form of quan- 
tum matter enabling in situ studies of most disparate 
physical phenomena, such as Josephson oscillations Q, 
or Anderson localization [1, Q , on a micrometer scale. 

In the limit of zero temperature and weak interac- 
tions, where all atoms of the gas can be assumed to share 
the same single particle quantum state, the dynamics of 
the condensate is described accurately by the nonlinear 
Gross-Pitaevskii equation . Finite temperature effects 
at thermal equilibrium are accounted for within higher 
order perturbation theories @, In contrast, only few 
theoretical works have been developed to model the non- 
equilibrium process of condensate formation itself. Pio- 
neering works, such as of Refs. [9l-[l2j. pointed primarily 
on the different dynamical stages of condensate forma- 
tion in terms of kinetic growth equations, and numerous 
efforts have led to highly accurate predictions for the time 
scale of condensate formation. 

Less is known about the condensate particle number 
distribution's dynamics in a dilute, weakly interacting 
Bose gas consisting of a fixed number of N particles. An- 
other question under current study fl3j is whether the 
equilibrium steady state of a dilute Bosc-Einstcin con- 
densate - which is finally reached only due to interatomic 
collisions, even in the case of very weak interactions - is 
unique, and characterized by thcrmodymanics and statis- 
tics of an ideal gas? And, how do quantum effects, such 
as number and energy fluctuations of the condensate and 
the non-condensate, which should become important es- 
pecially for mesoscopic Bose-Einstein condensates, evolve 
during Bose-Einstein condensation and eventually drive 
the Bosc gas into the final Gibbs-Boltzmann equilibrium 



state? 

Here, wc present a quantum master equation theory 
for a dilute Bosc-Einstcin condensate consisting of a fixed 
number of N particles. 

In contrast to previously derived effective equations for 
the condensate dynamics under the influence of the non- 
condensate environment, based for example on quantum 
kinetic theory (ll| or on an analogy with the laser master 
equation |14j , our condensate master equation fully takes 
into account conservation of the total number of parti- 
cles, i.e. of condensate plus non-condensate particles. In 
particular, the depletion of the non-condensate during 
the process of condensate formation results in conden- 
sate feeding and loss rates which arc different from the 
case where the Bose gas is coupled to an external particle 
reservoir with a fixed chemical potential. 

Apart from particle number conservation, our ap- 
proach relies on the separation of time scales between 
the condensate and non-condensate dynamics. The time 
scale To for condensate growth [L9l - |2l| is typically of the 
order of 1 — 4 s, and thus much slower than the timescale 
Tool ~ 10 — 50 ms of two-body collisions within the non- 
condensate [HI, [It], EH . We assume that these collisions 
lead to an effective thermalization of the non-condensate 
within each subspace of fixed non-condensate particle 
number, and to a rapid decay of non-condensate corre- 
lation functions with a rate T ~ r~j". Under these con- 
ditions, particle exchange between condensate and non- 
condensate leads to a Markovian master equation for the 
condensate's particle number distribution. Finally, we 
will show that its equilibrium steady state is given by a 
Gibbs-Boltzmann thermal state of non-interacting parti- 
cles under the condition that the Bose gas is sufficiently 
dilute, and that the decay of non-condensate correlations 
does not occur too fast, i.e. hf3T <C 1, where /3 = 1/fceT 
(fee being the Boltzmann constant and T the tempera- 
ture of the gas). 

The paper is organized as follows: The derivation of 
the quantum master equation is given in section[II] First, 
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wc define the condensate and non-condensate Hamiltoni- 
ans, and decompose all two body interaction terms in a 
physically motivated way. Using the assumptions men- 
tioned above, we then derive the quantum master equa- 
tion of Lindblad type for the reduced condensate density 
matrix. For dilute atomic gases in three-dimensional har- 
monic trapping potentials, the Lindblad master equation 
reduces to a simple rate equation for the condensate num- 
ber distribution, describing transitions No — >■ No ± 1 of 
the condensate particle number N . 

In Sec. IIII1 the rate equation is used to study the 
dynamics of condensate formation, and its equilibrium 
steady state. Wc numerically monitor the condensate 
particle number distribution during Bose-Einstein con- 
densation and extract time scales for condensate forma- 
tion. The steady state solution of the rate equation fi- 
nally yields a unique equilibrium steady state obeying de- 
tailed balance particle flow between condensate and non- 
condensate. In the case of weak interactions and under 
the condition Kf3Y <C 1 of not too rapidly decaying non- 
condensate correlation functions, the steady state turns 
into a Gibbs-Boltzmann thermal state of a canonical en- 
semble of N indistinguishable, non-interacting bosonic 
particles. 

We conclude in section HVl 



tion, 



II. QUANTUM MASTER EQUATION THEORY 

The separation of time scales between non-condensate 
thermalization and condensate growth motivates a de- 
composition of the gas into a condensate "system" and 
a non-condensate "environment" part, see Sec. Ill A[ In 
Sec. IIIB1 we then examine the two particle interactions 
between these subsystems as described by the Hamilto- 
nian, Eq. They fall into three different classes which 
we denote as single particle, pair and scattering events. 
Under the inclusion of all these two body interaction pro- 
cesses, the quantum master equation of Bose-Einstein 
condensation in a Bose gas with conserved particle num- 
ber N is finally derived in Sees. Ill Cl and ILTTJ1 
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which, as discussed in Ref. [8|, |22|, l23(, gives a good ap- 
proximation to the exact condensate mode at sufficiently 
low final temperatures, and sufficiently dilute atomic 
gases. 

In our following treatment, we will use "JoC?) as de- 
fined by Eq. (JT|) to describe the condensate wave func- 
tion also in a situation where initially not all particles 
occupy the condensate, and hence the condensate parti- 
cle number will change as a function of time. Neglect- 
ing the associated time dependences of ^o(r) is justified 
because it changes over the characteristic time scale for 
condensation, much longer than the time scale t co i. Wc 
can thus employ an adiabatic approximation and com- 
pute the rates in the master equation for a fixed 'J'o and 
at the end of the calculation only take into account their 
dependence on ^o-i hence on time. In the limit of very 
weak interactions, where the condensate state is approx- 
imated, at all times, by the ground state of the external 
trapping potential, |xo) ; or, if the initial state of the gas 
is already close to its equilibrium value, the situation is 
even simpler as the dependence can be entirely for- 
gotten. 

The total bosonic field expressed in an orthonormal 
basis {\^k),k G No}, where l^o) is the Gross-Pitaevskii 
ket, separates into 



* = 1*0)00 



E 



(2) 



with creation and annihilation operators Ofc and re- 
spectively, satisfying usual bosonic commutation rela- 
te, and [a k ,ai] = a\,a\ = 0. 



tions a,k , a\ 



2. Fock-Hilbert space 



A. Condensate and non-condensate subsystem 

After defining the condensate mode, we examine the 
decomposition of the full two body Hamiltonian in 
Eq. ([5]) into a condensate and a non-condensate part, 
and the interactions between them. 



1. Separation of the second- quantized field 

Quantitatively, we determine the condensate wave 
function 'l'o(r) (assuming all TV particles occupying the 
same condensate mode) by the Gross-Pitaevskii equa- 



The corresponding Fock states, forming a complete ba- 
sis of the many particle Hilbert space T on which these 
operators act on are denoted by \Nq, {iVfe}). The inter- 
pretation of a many particle Fock state is hence to find 
No particles in the condensate mode \^o), and {Nj.} = 
{Ni,N 2 , . . .} particles in the modes {|*i), |* 2 ), ■ ■ ■}• The 
basis {\^k),k <G N} is chosen such as to diagonalize the 
non-condensate Hamiltonian, see Eq. 1(5)). Let us point 
out briefly the tensor structure of the total Fock-Hilbert 
space T, corresponding to the subsystems condensate 
and non-condensate, respectively: 



T = To <& T± 



(3) 



As the condensate Hilbert space To is defined by To = 
span{|Ao) : No G N}, so is the Hilbert space T± of the 
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non-condensate by J-± = span{|iVi, N%, ...) : Nk G N}. 
Partial traces will be taken according to Eq. ^ in the 
following. 



3. Decomposition of the Hamiltonian 



The following decomposition of the Hamiltonian only 
requires the validity of the Gross-Pitaevskii equation for 
the condensate mode, and the orthogonality of the two 
fields an( i ^-L) m t ne sense that 



dr *J(r)#_L(r) =0 



(4) 



non-condensate fields and ^ T , , respectively: 



H± = J dr * ± (r) 



h 2 v~ 

2m 



+ V ext (r) 



*x(f) 



(8) 



where efe arc single particle energies of non-condensate 
particles. To model interactions between non-condensate 
particles, we assume that these lead to a rapid thermal- 
ization in the non-condensate thermal vapor, as will be 
further discussed in Sec. HID II below. 

Finally, the last term in Eq. ((6]), Vox, describes all 
two body interactions between condensate and non- 
condensate. This term will be examined in the following 
subsection. 



The Hamiltonian % in second quantization, including two 
body interactions Q, is given by 



U=J dr tft(r) 



h 2 v 

2m 



+ V cxt (r) 



(5) 



dr * t (r)* t (r)*(r)*(r) 



where \&(r) = \&o(*\) + *x(r) denotes the second- 
quantized bosonic field, and with g = Airati 2 jm quantify- 
ing the interaction strength in terms of the s-wave scat- 
tering length a. The neglect of three-body collisions im- 
plied by Eq. (0 is justified in the dilute regime ag 1 ^ 3 <C 1, 
with g the density of the atomic gas. The field decom- 
position in Eq. ([2]) splits the Hamiltonian H into three 
basic contributions, 



1~L — 7~Lo + 



(6) 



where Hq and H± describe a pure condensate and non- 
condensate, respectively. 

The condensate Hamiltonian Hq contains the single 
particle contribution linear in the field as well as the 
nonlinear, self-interacting two body term, and is given by 



U = Jdf $+(r) 



h 2 V~ 

2m 



+ V cxt (r) 



*o(?) 



dr * (r)* (r)*o(r)# (f) 



(7) 



When the average number of particles in the condensate 
is much larger than unity, a mean field approximation 
can be used to compute the ground state of Hq, which 
allows to recover the ordinary Gross-Pitaevskii equation 

©■ 

Concerning the Hamiltonian of the background gas, 
H±, wc first write down the contribution bilinear in the 



B. Two-body interaction processes 

Inserting the decomposition of the field ^(r), Eq. ©, 
into the Hamiltionian H, Eq. (O, we find, be- 
sides the condensate and non-condensate Hamiltonians, 
Eqs. ([3 [8]), various terms describing two particle interac- 
tion processes. Sorting these according to the number of 
condensate and non-condensate particles, which are cre- 
ated or annihilated during a two body collision event, we 
obtain 



where 



V~ 



Vxo = V_ + V_ + V D 



dr * 1 L (f)*^(r)*x(r)^ (F) 



g I dr *J(?)#^(rO#x (?)*_!_(?) 



(9) 



(10) 



accounts for single particle events, where the conden- 
sate particle number changes by AiVo = ±1, and cor- 
respondingly, the number of non-condensate particles by 
AiV ± = Tl- 
Second, 



| / dr *5_(r)*J_(?)*o(r)*o(iO 



,t/* 



(11) 



describes pair events, where two condensate particles are 
created or annililated, i.e., AiVo = ±2 and AN± = =p2. 
Finally, the term 



V?j = 2g dr ^(r)* (r)* ± (r)* (r) 



(12) 



describes scattering events, where the number of conden- 
sate and non-condensate particles is unchanged (AiVo = 



4 




k* 



I* 



/ \ 

o o 




ficiently low temperatures (i.e., a sufficiently peaked con- 
densate number distribution close to TV). Indeed, when 
we combine upper left diagrams in Fig. [TJ and their her- 
mitian conjugates, with mixed single particle contribu- 
tions in Eq. ([5]), we get the vanishing term 



dr * T 



2m 



+ y cxt (r) + #T(r)* (r) 



*o(f) 
(13) 



In total, the Hamiltonian T~L in Eq. ([5]) thus decomposes 
into 



T~L — Wo ~\~ Wj 



V 



(14) 



FIG. 1: Diagrammatic representation of all two body loss pro- 
cesses m Eq. ©. The upper two diagrams show single particle 
losses (~*), where one non-condensed particle is created and 
one condensate atom is annihilated. The upper left diagram 
denotes a term linear in the non-condensed field, 0('&±), 
which vanishes in combination with crossed single particle 
terms as a consequence of the Gross-Pitaevskii Eq. (JTJ and 
the orthogonality condition in Eq. Q. The lower diagrams 
show pair losses («™», lower left) and scattering processes (O, 
lower right). The conjugate processes (not shown), related 
to condensate feeding, are obtained by exchanging the corre- 
sponding labels with respect to the diagram center. 



ATVl = 0). As we will see later, scattering events do 
not contribute to the master equation for the condensate 
density matrix, which will mainly be governed by single 
particle events, with negligible influence of pair events. 

To illustrate the different interaction terms, we intro- 
duce a diagrammatic representation of the interaction 
matrix elements [IB, Hit- These are depicted in Fig. [T] 
Annihilation and creation of condensate particles are 
denoted by and 0*, respectively, whereas k,l,m and 
k*,l*, m* refer to annihilated, or created particles of the 
corresponding non-condensate modes. Note that Fig. Q] 
contains only condensate loss events, where the number 
of condensate particles decreases. The conjugate pro- 
cesses, corresponding to condensate feeding, arc obtained 
by exchanging the corresponding labels with respect to 
the diagram center. 

Furthermore, Fig. Q] also shows processes of first order 
in the non-condensate field (upper left diagram), which 
are, however, not contained in Eq. ([9]). The reason is that 
these processes cancel out with mixed, single particle con- 
tributions between condensate and non-condensate fields 
in the Hamiltonian in Eq. ([5]). This is a consequence 
of the orthogonality of the two fields, *J(r) and 4 , j_(r), 
see Eq. (j4]), and the fact that v I / o(f) is an approximate 
solution of the Gross-Pitaevskii equation, Eq. ([T]), for suf- 



where the different interaction terms in Eq. (JTOJ) , 

in Eq. (fTTj) , and Vq in Eq. (|T2"|) , account for single particle 

(~*»), pair (<~^>) and scattering (O) contributions. 



C. Evolution equation for the total density matrix 

In analogy to the standard quantum optical deriva- 
tion [Tol llq ]. we start with the von- Neumann equa- 
tion, considering a many particle state <xW (t) of fixed 
particle number TV, defined on the Fock-Hilbert space 
J r = J 7 ®J 7 ± in Eq. (©: 



dt 



(15) 



where T~L is the total Hamiltonian, see Eq. With 
the decomposition of H. in Eq. ((6]), the von-Neumann 
equation turns into 



dt 



H ,<J iN) (t) 



n u o {N \t) 



(16) 



Note that we use here the linearized non-condensate 
Hamiltonian in Eq. ([8]). We transform all operators, i.e., 
the condensate and the non-condensate field, ^o(r) and 
^j_(r), as well as the density matrix a^- N \t), to the inter- 
action picture (denoted by the label J), which is carried 
out with respect to the Hamiltonian parts Hq and TL± 
of the subsystems condensate and non-condensate. The 
different operators hence undergo the transformation 



X(t) ->JfW(t) =U{t)XU\t) 



(17) 



with respect to the time evolution operator U (t) given by 



U{t) = exp 



-(H +Hx)t 



(18) 
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The time evolution of the full density operator cr^ 7 ) (£) 
in the interaction picture is then determined by the inter- 
action between condensate and non-condensate particles, 
according to: 



BOSE GAS 



da^^jt) 
dt 



* W *>(t)] . (19) 



where V^{t) is obtained by inserting the time depen- 
dent annihilation (and creation) operators, e.g. ao(t) = 
aoe~ ltJ/ot / h and ak{t) = a,ke~ ltkt l h , in the corresponding 
time independent expressions derived in the previous sec- 
tion. Integration of Eq. ([T9"]) between t and t + At leads 
to 



&^ T \t + At) 









t+At 


i 


/di 


h 






t 



(0 



mt'u^if) 



(20) 



For short times, At, a good approximate solution of 
Eq. (|20[) is obtained by its iteration up to second order in 
Vq± (which is required since the first-order terms vanish, 
as we will see later): 



t+At 



Aa^(t) 



i 



dt' 



v^iV), 



(t) 



t+At t' 

> / r dt " 



(21) 



where we have set Aa^ N ^{t) = a^ N ^(t + At)-^ N ^(t). 
Note that Eq. (|2"Tj) expresses the state at time t + At 
(left-hand side) fully, as a function of the state a^' 1 ^ (t) 
at time t - in contrast to the exact Eq. (f20f , where states 
( j( Ar ^ 7 )(i') at all intermediate times t' appear on the right- 
hand side. 



D. Time evolution of the reduced condensate 
density matrix 



The time evolution of the condensate in the presence of 
the non-condensate gas is obtained by taking the partial 
trace over F±_ in Eq. (j2"Tj) . To get a Markovian mas- 
ter equation for the reduced condensate density matrix, 
p (t) = Trj^ ± a^ N \t), we use a Born-Markov ansatz 
generalized for the TV-particle state a^ N '(t) which allows 
to express a^ N '(t) completely in terms of the reduced 
condensate density matrix p Q N \t) at time t, see Eq. (f22j). 



1. N on- condensate thermalization 



In standard derivations of master equations for systems 
coupled to thermal reservoirs 

[H El, the Markov as- 
sumption is justified by assuming a thermal state pe(T) 



energy 



•x 

particles 













EXTERNAL 
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RESERVOIR 











condensate 



non-condensate 



FIG. 2: Schematics of microsopic many particle dynamics. 
The total number of atoms in the Bose gas is fixed to N 
and conserved during condensate formation. Atomic colli- 
sions within the non-condensate are modelled by coupling the 
non-condensate part of the gas to a heat reservoir which is at 
fixed temperature T. The condensate part is initially not at 
equilibrium with the non-condensed fraction, and both sys- 
tems undergo a net exchange of particles, induced by atomic 
two body collisions between condensate and non-condensate 
atoms, which are fully taken into account in the derivation of 
the master equation. The final equilibrium steady state of the 
gas exhibits detailed balance particle flow between condensate 
and non-condensate. 



for the environment, which is supposed to be unchanged 
by the interaction with the system. Then, the total state 
a(t) would be given by the product &(t) = po(t)® Pe(T), 
hence completely determined by the reduced state po{t) 
of the (condensate) subsystem. However, in our case, 
this simple product ansatz cannot be applied, since con- 
densate and non-condensate are correlated by particle 
number conservation: If one finds Nq particles in the 
condensate, the particle number in the non-condensate 
is determined as N± = N — No, and vice versa. 

The physical origin of the non-condensate thermaliza- 
tion is the interaction between non-condensate particles, 
which leaves the number of non-condensate particles un- 
changed. We hence couple the non-condensate to a heat 
bath only allowing for exchange of energy, but not of par- 
ticles, see Fig. [5] The thermalization then occurs only 
within subspaces of fixed particle number. In addition, 
we assume that coherences between subspaces of differ- 
ent particle number are destroyed due to the coupling 
with the heat bath. Under this assumption, the total 
iV-particle state is obtained as: 

N 

&W(t)= J2 PN(N a ,t)\N )(N \® P± (N - N ,T) , 

2V =0 

(22) 

where pN(No,t) denotes the probability of finding Nq 
particles in the condensate (or, equivalently, N — Nq par- 
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tides in the non-condensate), and 

s /at » r Q.N-N e~ l3 ' H± Qn-n ,„„x 

PJ .(^-JVb,T)= Zx{n _ NqiT) (23) 

describes a thermal state projected onto the subspace 
of N — Nq non-condensate particles, with corresponding 
projector Qn-n i an( i normalization factor [2a ] 



2j.(JV-JV 0j r) = 'Erjr L {q 



Qjv- 



iv } • (24) 



Note that, since (t) is diagonal in the Fock basis, it 
is invariant under the free evolution U(t), Eq. (fT5)l . and 
hence a^ N,I \t) = a^ N ^(t). In the following, we hence 
drop the index 'J' referring to the interaction picture for 
the iV-particle state a^ N '(t), or its reduced condensate 
state p Q N \t), sec below. 



2. Evolution equation for the condensate density matrix 

Taking the partial trace over the non-condensate, we 
obtain the reduced condensate density matrix: 

N 

p { N) (t) = Trr x {.(">(()} = Yl PN(No,t)\N )(N \ . 



JVn=0 



Obviously, also the reduced condensate state is diagonal 
in particle number representation as a direct consequence 
of our assumptions on particle number conservation and 
rapid non-condensate thermalization. Thus, both, the 
reduced condensate density matrix, Eq. (|25p . as well as 
the total iV-particle state, Eq. (|22|) . are completely de- 
termined by the condensate particle number distribution 
PN(N ,t). 



(25) 



Inserting Eq. ([22]) in the right-hand side of Eq. (|21j) . 
and taking the partial trace over the non-condensate, 
leads to a closed evolution equation for the reduced con- 
densate density matrix. Moreover, it can be shown that 
the terms of first order in the interaction V<^ vanish after 
taking the partial trace over Indeed, from the diag- 
onal form of the iV-particle state a^(t), see Eq. (|22|). it 
follows that 

Tr^{[v^(O^ (ArJ) (i)]}=0. (26) 
From the remaining second order terms, we obtain: 



N 

£ 

N =0 



t+At 



Vox (*') . ^ox (*") . Pn (No ,t)\N Q )(N \ (8 p± (N - N , T) 



(27) 



Writing the interaction term (t') as a sum over the 
three different processes (single particle, pair and scat- 
tering events), see Eq. ([9]), we can now verify that any 
mixed commutator in Eq. (|27[) vanishes - again as a con- 
sequence of the diagonahty of a^ N '(t). 

Hence, single particle, pair and scattering events in 
the gas are dynamically independent from each other. 
Furthermore, it can be shown that scattering events, de- 
scribed by Vq, do not contribute, since they leave the 
number of condensate particles unchanged. We are left 
with: 



A/T(t) _ Ap^W 



At 



At 



A/T(*) 



At 



(28) 



where the two terms on the right-hand side of Eq. 



are obtained by inserting the corresponding interaction 
terms VL» and VV~>, instead of the full interaction V^o 
into Eq. ([27)1 . 

3. Quantum master equation of Lindblad type 



In order to perform the time integration in Eq. (|27l) . 
we first notice that the right-hand side depends only on 
the time difference r = t' — t". Second, we assume that 
only times r <C At contribute to the integral, due to 
the rapid decay of non-condenate correlation functions. 
To implement this rapid decay, we assume that the two 
point correlation functions of the non-condensate decay 
on the average time scale t co \ of a two body collision 
event. Performing the time integral as 
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*t+At pt' pAt pt+At pAt 

/ dt dt" = dr / dt' -At dr , 

Jt Jt JO Jt+T JO 



(29) 



using t < At for the second equality. However, even 
though At has to be larger than the decay time ~ 
r co i of non-condensate correlation functions (see below) , 
it is still much smaller than the time scale Tq for the 



r 



condensate evolution. In this case, the coarse-grained 
rate Apo(t)/At can be replaced by the instantaneous 
time derivative dpo(t)/dt to obtain the following Lind- 
blad master equation: 



^OT 1 - E ^(iV ,T) [4(iVo) / 5^(t)«S+(iV ) - i {<Sj(JV )4 (^o), (*)} 

JV = 0, L + 

3=+.- 

N r 

5] 7^o,T) V j (N )pi N \t)v}(N )--{p}(N )V j (N )ji N \t)} 



N =0, 
3=+,- 

where the quantum jump operators <S±(iVo), and T > ±(Nq) are defined by 

S+(N Q ) = \N Q + 1)(N \, S_(N ) = \N a - 1)(N \ 
V+(N ) = |JV + 2><JVo|, V-(N ) = \N Q -2)(N \ . 

I 



(30) 



(31) 



Obviously, S+(Nq) adds one particle to the condensate 
with a rate $(N ,T) = 2{N + 1)A+ (N - N ,T), 
whereas S_(N ) destroys a condensate particle with the 
rate Q(N ,T) = 2N Xz t {N -N ,T), given a number of 
(N — No) non-condensate particles, and a temperature T 
of the heat reservoir. In a similar way, V±{Nq) describe 
the simultaneous creation of two condensate particles 



with a rate 7 + (N ,T) = (N + l)(N +2)\+ t (N-N , T), 
and the depletion of two condensate particles with a rate 
J N (N ,T) = N {N - l)A_(iV - Nq,T). The different 
transition rates ^(Nq,T) and ^(No,T) are defined by 
the following integrals over non-condensate correlation 
functions: 



|J ^ dr df' ^(f}* (r) J dr e^V^^F, F, N - N , T, r) 

( OO 

\t>(N-N ,T) = Re| //d?df * (f)*o(f)*S(f )*o(0 / dr e^V^^F, F, iV - 7V ,T,r) 



>33) 



where F\iV - Nq,T,t) and r,iV - 

Nq,T,t) are correlation functions of the non-condensate 
field for single particle (~-») and pair (+~*) events, given 
that (N — No) particle are in the non-condensate gas. In 
Eqs. (|32l [33]) . we have extended the time integral from 
At to oo, assuming a Gaussian decay of non-condensate 



correlations due to thermalization which occurs within a 
time interval on the order of the average time T = t^ 1 
for two-body collisions. 

The remaining coherent parts Q^' and of non- 
condensate correlations are determined by the thermal- 
ized state in Eq. 
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g(+)(r, r', N - N , T, r) = r)*l(f, r)<F ± (r, r)*^(f', 0)* ± (f', 0)* ± (f', 0))^ ^ 

Q# (r, r', JV - N , T,t) = { r)* ± (?, r)¥|_(?, r)^(/, O^f, 0)tf ± (f, 0) ) 

for single particle events, and by 

^(r,f,iV-iVo,r,T) = (^(r,r)^(r,r)* ± (f', 0)<F ± (r', 0) ) 

. . / « . « „ . v (N-N ) 

g£(?,?,N-N 0) T,T) = (*x(r,r)^(?,r)*l(r',0)*i(f,0))^ 



(34) 



(35) 



for pair events. In Eqs. (|34l [35)) . (. . .)nv_jv ) denotes 
the average Tr^ {. . . p±(N — No, T)} with respect to a 
thermal non-condensate state with (N — No) particles. 
Note that the imaginary parts of Eqs. (|32H33p which 
lead, in principle, to a shift of the condensate energy 
levels (similar to the Lamb shift known from quantum 
electrodynamics [HI]), drop out from the master equation 
due to the diagonal form of the reduced density matrix, 
Eq. ((251). 



4- Quantum master equation of Bose-Einstein condensation 

From the master equation of Lindblad type in 
Eq. (|3T)|) . we can derive the evolution equation for the 
condensate particle number distribution, pjv(iVo,£) = 
(No\p N \t)\No) . Considering only single particle pro- 
cesses (~~>), since they dominate the condensation process 
in three-dimensional harmonic traps, see section III D 51 
leads to the quantum master equation for quantum jump 
processes with ANq = ±1: 



dp N (N ,t) 
dt 



[t+(N ,T) + t-(N ,T)}p N (N ,t) 



+ Z+(N -l,T) PN (N -l,t) 
+ £ N (N + l,T)p N (No + l,t) 



feeding losses 

|No+lxNo+l| 



S(+.N<.) 



^(-N„+l) 




£(+,No-l) 



S(-.No) 



(36) 



FIG. 3: Probability flow between different condensate num- 
ber states as expressed by the quantum master equation 
in Eq. (|36[) . The corresponding transition rates are de- 
fined by &(N ,T) = 2(N + 1)X+(N ,T), and ^ N (N ,T) = 
2N \~{N ,T), with \t(N ,T) given by Eq. In 
the stationary state which is reached for long times t — s> 
oo, the rates obey the condition of detailed balance: 
Z+(N , T)p N (N , T) = Q(N + 1, T) PN (N + 1, T). 



the steady state of the system is therefore reached, if, and 
only if the net probability flux for every state |JVo)(JVo| 
(— >• detailed balance of probability flow) is zero, i.e. 
Z+(N ,T)p N (N ,T) = Q(N + l,T)p N (N + 1,T) for 
aU7V = 0...N. 



with t+(N Q ,T) = 2(N + l)\+(N - N ,T), and 
£jf(N ,T) = 2N \Z t {N - N ,T), where the transition 
rates A± (N - N ,T) are given by Eq. (J32J. 

Bosc-Einstcin condensation is now reduced to a simple 
rate equation, the master Eq. (|36p . which describes in 
particular the buildup of a macroscopic condensate occu- 
pation from the fluctuating thermal vapor. As sketched 
in Fig. [31 net particle flow towards a state |iVo)(iVo| is 
described by the current £n(Ao — 1,T)pn(Nq — l,t) + 
£n(N + 1,T)pn(No + l,t), and particle flow from the 
state \N )(N \ by the current £+ (N , T)p N (N , t) + 
£x(N ,T)p N (N ,t). As will be shown in section IfflDl 



5. Transition rates for Lindblad dynamics 

We now evaluate the different decay rates, Eqs. (J35J- 
|33|) . For this purpose, we decompose the higher order 
correlation functions of the non-condensate fields into 
second order correlation functions according to the Wick 
theorem [26|, [27| , and perform the integrals over r , r' and 
r, see appendix [XJ For the single particle creation and 
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loss events, the result is: 
16Tr 3 h 2 a 2 



[f?(k,l,m,N ± ,T)x 



X (5 (r) (tJfc + L0 t - UJ m - Wo) + 

+g^(k,l,m,N ±1 T)S^(LJ l -Lj ) 

where uj k = e k /H, u>o = and where 

f?(k,l,m,N x ,T) = (iV fc )(iV ± ,r)(7V i )(7V ± ,T)x 



x [(JV TO )(JV X> T) + 1]|C 



and 



g+(k, I, m, Nj_, T) = 2(N k )(Nj_, T) (N) (N^T)x 
x (N m )(N ± ,T)C k fC^ 



(37) 



(38) 



(39) 



are the weight functions for condensate particle feedings. 
Correspondingly, 



and 



/r (A, I, m, iV±, T) = [(JVfc)(iVx,T) + 1] x 
x [{N,)(N ± ,T) + 1] (N m )(N ± ,T)\Q°\ 2 

gT(k,l,m,N ± ,T) = 2(N k )(N ± ,T)x 

x [(JV,)(M,T) + 1] (N n )(N ± ,T)CCot 



(40) 



(41) 



are the weight functions for condensate particle losses. 
The functions f±(k,l,m, N±_,T) and g±(k,l,m,N±,T) 
depend on temperature, the number of non-condensate 
particles, N±_ = (N — No), and on the quantum mechan- 
ical probability amplitudes 



/•Ik, /AmOu 



d?*S(?)*™(?)^(?)*/(?) (42) 



for the different microscopic single particle feeding and 
loss processes with energy balances Aw^ = Ae^/h = 
(u) k + u>i — LJ m — ujo) or (lui — wo). The ^-distribution 
in Eq. (37J), S^(A^) = 0F/r 2 cxp[-(A^) 2 /4r 2 ], re- 
flects conservation of energy during the different single 
particle feeding and loss processes on a certain width \/2r 
arising from the decay of the non-condensate field corre- 
lation functions in Eqs. (|3"2l 13"3"1) . Therefore, only single 
particle processes with energy balances Ae^/h < T will 
contribute to the rates in Eq. (|37l) . 

The average occupation number (N k )(N± , T) of a non- 
condensate single particle state \^k), given a thermal 
state projected onto the subspace of N±_ = (N — No) 
non-condensate atoms, reads (see appendix IB"]) 



(N k )(N ± ,T) 



(43) 



exp[/?(e fe -/xx(iVx,T)]-l ' 
where (N± , T) is defined by the normalization con 



dition: J2 k¥l0 (N k )(N±,T) = N±. According to this 
definition, fx±(N±,T) equals the chemical potential of 
a thermal state of N± = (N — No) non-condensate parti- 
cles [ID]. From (N k ) + 1 = (N k ) exp[/3(e/-— /ax)], and using 
the energy conservation as expressed by the (S-function in 
Eq. (|37[) . one can derive the following relation between 
the single particle loss and feeding rates: 

A+(iV - N ,T) = e^^-^'^A" (iV - N ,T) , (44) 

where A(i(N - N ) = fi±(N - N ,T) - /Lt . To obtain 
Eq. (|4"4"| . the finite width V of the (5-function is neglected, 
which is justified under the condition HTf3 <C 1. The 
relation (|44|) will be useful in Sec. IIII fJl to determine the 
equilibrium state of the Bose gas. 



Note that Eq. ([4"3"| takes into account the depletion 
of the non-condensate during condensate formation, en- 
suring that (Nk) — >• as iVo —> N. According to 
Eqs. (|38I40[) . also the condensate feeding and loss rates 
tend to zero in this limit. In contrast, the rates obtained 
within quantum kinetic theory [TTj increase with increas- 
ing condensate particle number No- 



Finally, the rates for pair events turn into: 



J2 f?(k,l,N x ,T) 



x 5 (r) K +ui -2wq) , 



(45) 



with the weight function 

fT(k,l,N ± ,T) = (N k )(N ± ,T)(Ni)(N ± ,T)\ffi\ 2 (46) 
for pair feedings, and correspondingly 

f-(k,l,N ± ,T) = ((N k )(N x ,T) + l)x 



x((JV,)(iVx,T) + l)|&°| 2 
for pair losses, with 



^00 " 



dr (*g(r))** fc (r)*,(f) . 



(47) 



(48) 



Looking at the energy balance of a pair event, Ae^/h = 
uj k + lui — 2wo, we see that pair events occur as energy 
non-conserving processes, i.e., Ae^/h^> T ~ 10— 20 Hz, 
since the single particle energy fiQ of condensate particles 
is smaller than the energies tk.i 3> hT of non-condensate 
particles in a thrcc-dimcnsonal harmonic trap. There- 
fore, pair events can be neglected in comparison with the 
single-particle events in the master equation (|36|) . For 
the same reason, we can neglect the terms associated to 
the g± functions in Eq. (|3T|) as compared to those asso- 
ciated to the f± functions. 
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FIG. 4: (color online) Time evolution of pN{No,t) (probabil- 
ity is indicated by the color gradient in the figure) according 
to Eq. (|36[) during condensate formation for parameters sim- 
ilar to Ref. 0], i.e. in a gas of N = 2000 87 Rb atoms con- 
fined a harmonic trap with frequencies u x = ui y = 2n x 42.0 
Hz and u) z = 2n x 120.0 Hz. The temperature has been 
set to T = 20.31 nK, and the critical temperature of the 
gas is T c — 33.86 nK. According to the atomic density g — 
2.6 x 10 12 cm -3 and the s-wave scattering l ength a = 5.7 nm, 
we set F = (a 2 g)~ 1 v = 34 Hz, where v — ^/3fcsT /m. 



III. RESULTS 

In this section, we present numerical studies of the 
condensate particle number distribution obtained from 
Eq. (|36|) during Bose-Einstein condensation, and derive 
the unique equilibrium steady state of the Bose-Einstein 
condensate. Under the assumption hT(3 <C 1, the equi- 
librium steady state is proven to be a Gibbs-Boltzmann 
(thermal) state of non-interacting particles in the dilute 
regime aN/L <C 1, with L = h/muj the extension of the 
harmonic oscillator ground state. 

A. Perturbative calculation of transition rates 

For this purpose, we now consider the case of very 
dilute, weakly interacting gases. Since the transition 
rates derived in the previous section originate from pro- 
cesses of second order in the interaction Vj_ , all the rates 
are proportional to a 2 , as evident from the prefactors in 
Eqs. (|3"TI B5T) . The remaining dependence of the rates on 
the interaction strength originates from the single par- 
ticle wave functions \^>k), the eigenvalue of the Gross- 
Pitaevskii equation an d the non-condensate single 
particle energies which are themselves functions on 
the parameter ag, see Eq. ((T|). 

Interested in the case of dilute and weakly interacting 
gases, we thus expand these quantities (i.e. I^o)) A*o and 
in terms of the scattering length a, taking into ac- 
count only the first non- vanishing contribution, given by 
their non-interacting limits. Hence, the basis states \^k) 
turn into the single particle eigenstates \\k) of the trap- 
ping potential, with corresponding ground state energy 



eo, whereas the e' k s are the energies of the excited states. 
Thereby, we evaluate the transition rates to lowest non- 
vanishing order - proportional to a 2 - in the s-wave scat- 
tering length. 

Quantitatively, this procedure is correct as long as the 
ground state of the Gross-Pitaevski equation ([1]) is well 
approximated by the single-particle ground state. This, 
in turn, is the case if the interaction energy g./VI'J'ol 2 is 
much smaller than the harmonic oscillator energy hui, or, 
in other words, if aN/L <C 1, where L = H/muj denotes 
the extension of the harmonic oscillator ground state. 

When aN/L increases, the ground state of the Gross- 
Pitaevskii equation is progressively distorted, as well as 
the various excited states, making the explicit calcu- 
lation of the rates, Eqs. (|37H42]) , more difficult. How- 
ever, no drastic change is likely to take place, making 
the following predictions qualitatively and maybe semi- 
quantitatively correct for a realistic situation like Bosc- 
Einstcin condensation of a Rb or Na gas. 



B. Dynamics of Bose-Einstein condensation 

Equation (|36p is solved numerically to propagate the 
condensate particle number distribution pN(No,t) in 
time. Figure [5] displays a typical example for the time 
evolution of pn{No, t) for a gas of N = 2000 87 Rb atoms 
which undergoes the Bose-Einstein condensation phase 
transition in a three-dimensional harmonic trap with fre- 
quencies uj x = u) y = 2n x 42.0 Hz, w z = 2tt x 120.0 Hz. 
The final temperature of the gas is T — 20.31 nK, given 
an ideal gas critical temperature T c = 33.86 nK [||. Note 
that, r is a free parameter in our theory, provided it is 
smaller than ksT and larger than the external trap fre- 
quency. However, we have numerically checked that the 
transition rates do not significantly change with varying 
r in this parameter regime. 

To calculate the transition rates leading to the conden- 
sate growth scenario in Fig. [U we used the semi-classical 
limit 0, [2l| , where the discrete sums in the feeding and 
loss rates in Eq. (|37|) are replaced by an integral over the 
density of states 17(e) = e 2 /2(h 3 uj x u>yOJ z ). This shifts the 
final condensate fraction by appr. 10% as compared to 
the exact numerical evaluation of the discrete sums (em- 
ployed in Figs. [5] and [J5] below), but does not change the 
qualitative behavior observed in Fig. 2J 



C. Average condensate growth 

From the time evolution of the distribution pN(N ,t) 
defined by Eq. ([3^]) the growth of the average condensate 
population can be extracted using 

jv 

(N Q )(t) = Y, N oPN (N ,t) , (49) 

JV o =0 
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FIG. 5: (color online) Comparison of condensate particle 
number distribution arising from the master equation (red 
solid line) vs. the Boltzmann thermal state ansatz of Eq. (|54p 
(blue squares), for the same parameters as in Fig. [4] The gas 
temperature is T = 25.0 nK, i.e. /3KT ~ 0.01, with T = 34 
Hz. 



and deriving a corresponding condensate growth equa- 
tion (Tl| . For this purpose, we assume a sufficiently nar- 
row peaked distribution p^(No, t) around the mean value 
(No) as indicated by Fig. 21 such that the rates are ap- 
proximately constant in this narrow region, meaning that 



Xt (N - No, T) a Xt (N - (N ), T) 



(50) 



for No close enough to (No). Taking the average 
N dp(N ,t)/dt with dp(N ,t)/dt given by Eq. (J36j) fi- 
nally leads to the following growth equation for the av- 
erage condensate occupation: 



d(N ) 
dt 



= t+({N ),T)-£j f {(N ) + l,T) 



(51) 



with £+((No),T) = 2((N ) + 1)A+ (JV - (N ),T), and 
Q((No) + 1,T) = 2«JVb> + 1)A- (JV - (JV ),T). The 
equilibrium value of (JVo) is hence defined by the detailed 
particle balance condition At(JV - (JV ),T) = A^(JV - 
(No),T). According to the above relation between the 
rates A+ and A^, sec Eq. f4"4"|) . this implies equality of 
the chemical potentials on average: 



Ho = H ± (N-{N ),T) 



(52) 



In the next section, we will show that not only the average 
condensate occupation, but also the whole steady state 
distribution pn( No, T) agrees with the thermodynamical 
prediction. 



D. Steady state distribution 



For this purpose, we solve Eq. (|22)) for the steady state 
distribution, defined by dp^(No,t)/dt = 0, which leads 
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FIG. 6: (color online) Average condensate fraction (No)(t)/N 
(upper panel), and standard deviation ANo (lower panel) of 
the condensate particle number distribution pjv (JVo , T) as a 
function of temperature, and otherwise the same parameters 
as in Fig. [5] The master equation's steady state distribution 
(red solid line) agrees very well with the canonical ensemble 
prediction (blue squares) of Eq. (|54|l . 



to: 



(53) 



Let us now compare this steady state to a thermal state 
of JV non-interacting particles at temperature T: 



Z(N, T) ' 



(54) 



with the partition function Z(N,T) of JV indistinguish- 
able particles, and the projector Qm onto the Fock space 
of JV particles. In the absence of interactions, T~L in 
Eq. (|5~4")) is the Hamiltonian of the gas in Eq. ([5]), with 
g = 0. To proof the equality of the state &N,th and the 
steady state of the Bose gas in Eq. ([22]) . with pn(Nq, T) 
given by Eq. (|5"5|) . it needs to be shown that the recursion 
relation for the condensate particle number distribution, 



PnmANq, T) 
PN,th(N + l,T) 



^„ Z ± (N-N ,T) 
Z x (N-No-l,T) ' 



(55) 



which results from tracing Eq. (|54p over the non- 
condensate, applies as well for the steady state, Eq. (|53|) . 
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of the master equation. In Eq. (|55j) . Z±(N — Nq,T) is 
the partition function of (N — N ) non-condensate parti- 
cles, see Eq. (|24|) . and eo the single particle ground state 
energy of a non-interacting gas. This can be seen if we 
approximate A~ (iV - N ,T) « A^(iV - N + 1,T), ne- 
glecting terms of the order of N~ 1 . In this case, we obtain 
from Eq. ((551) : 

Pn(N ,T) ^ X-(N-N ,T) = « (eo _^(iv X) T)) 
PN (N + 1,T) \+(N-N ,T) 

(56) 

where we used Eq. (j44[) . and eo = Mo ( m the regime of 
small interactions, aN/L <C 1, see Sec. IIII A[) . 

Now, the non-condensate chemical potential, as de- 
fined b y th e normalization condition in Eq. (|43|). can be 
shown [25J to be related to the non-condensate partition 
function, via -(3[i±(N - N ,T) = \nZ ± {N - N ,T) - 
lnZ^(7V — A*o — 1,T). Therewith, we arrive at the re- 
currence relation in Eq. (|55[) , which was to be shown. 
Hence, the steady state of the entire Bose gas in Eq. ([22)) 
is given by the thermal state in Eq. (|54")) , in the case of 
weak interactions, and under the condition /3hT <C 1 for 
which Eq. (|4"4")) is proven to be valid. 

The 1/iV approximation required for the above proof 
is confirmed by comparing the exact numerical calcula- 
tion of the steady state condensate particle number dis- 
tribution to the prediction of the Boltzmann ansatz in 
Eqs. (fMf . Fig. [5J shows the stationary particle number 
distribution for the same parameters as in Fig. 01 In or- 
der to show that the agreement holds up to the critical 
temperature (and beyond), Fig. [5] displays the compari- 
son of average condensate occupations, and the standard 
deviation of the stationary condensate particle number 
distributions (such as the one depicted in Fig. [5j, as 
a function of the entire range of relative temperatures, 
T/T c , for N = 2000 atoms for the same trap parameters 
as in Fig. [5j Again, we observe close agreement between 
master equation and the Boltzmann ansatz: The shift of 
the critical temperature is about 10% with respect to the 
critical temperature T c of Bosc-Einstcin statistics in the 
semiclassical limit m in both cases. 



IV. CONCLUSION 

We have presented the conceptual part and first nu- 
merical results of a number-conserving quantum master 
equation theory to describe the transition of a dilute gas 
of N bosonic atoms into a Bosc-Einstcin condensate. The 
central result of our theory is the quantum master equa- 



tion in Eq. (f36|) which describes the time evolution of 
the reduced condensate state in contact with the non- 
condensate environment for a fixed total atom number. 
In the dilute gas regime, we numerically monitored the 
full condensate particle number distribution piy(No,t) 
during condensate formation. 

The theory predicts condensate formation times of the 
order of seconds, matching experimentally and theoret- 
ically observed times scales [3, 0. The derived steady 
state for a dilute, weakly interacting Bosc-Einstein con- 
densate undergoing Markovian dynamics is unique, and 
proven to obey the same statistics as a Gibbs-Boltzmann 
thermal state of non-interacting particles, in the case of 
weak interactions aN/L <C 1, and for the case H(3T <C 1 
of not too rapidly decaying non-condensate correlation 
functions. 

Future improvements of our model will consist in a 
microscopic derivation of the rate T describing the de- 
cay of non-condensate field correlation functions (e.g. by 
diagrammatic expansion techniques for higher order cor- 
relation functions), which, in the present version, has 
been introduced in a rather phcnomcnological way. Fur- 
thermore, the condition aN/L C 1 of weak interactions 
may be relaxed. Since, in this case, the condensate wave 
function will depend on the number of condensate par- 
ticles, this will in particular require to introduce time- 
dependent condensate and non-condensate wave func- 
tions. Finally, it remains to be studied whether devia- 
tions from the Gibbs-Boltzmann occur if the condition 
h/3T < 1 is not fulfilled. 
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Appendix A: Two point correlation functions 

Here, we decompose the correlation functions of of the 
non-condensate field for single particle (~+) and pair (<-~->) 
processes into products of two-point correlation functions 
with Wick's theorem, which applies to thermal expecta- 
tion values [MUzl- 

We begin with the correlation function for single par- 
ticle processes </(+>(?, r,N- N , T, r) in Eq. JMJ): 
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g(+) (r, ?, N - N , T, t) = (¥ x (r, r)¥ ± (r, r)* ± (r, t)& x (?, 0)$x(r , 0)$±(r, 0) ; 

.(Mo)/ '-N ) 



(N-N ) 



2 r)* x (f,0) 



(JV-iVo) 



{N-No) , 

*l(?,r)tf L (?' 1 0)) (* ± (? ! r)*l(r',0) 



4(^(r,r)* ± (r,r)) V " U '(£ x (r, r)*x(r , 0)V"' (£ x (r, 0)*x(f, 0) 



(N-N ) 



Fx 



(N-Nq) 



(Al) 



and similarly: 



g^)(r,r',iV-iVo ! T,r) = ( £ x (r, r)* ± (r, t)* ± (F, t)&J? , 0)^(f, 0)#x(r, 0) ; 

. (JV-ATo) / „ , \ (N-N ) I „ „ . . \ (JV-iVo) 

ft ^iTf . rs' rv\\ /iTi , in ^\,T ( t rs / rrA /iiv . at ^ - 



(JV-JVo) 



2<*i(?,r)*x(? J 0) 



■Tx 



*i(?,t)*HF,0) 



*_l(F, T)* 1 L (f ,0) 



\ (N-N ) I „ » . \ (N-N ) I „ . 

4(#5_(r,T)*x(? ) r)) (*x(?,t)*5_(^,0)) (*+.(?', OJ^xC^O) 



(AT-iVo) 



(A2) 



The non-condensate field $i(f,T) in the interaction picture with respect to H± in Eq. l[8]). written in the single 
particle basis set e N}, turns into 



*x(r,r) = Wx(r)*x(r)Wl(r) = ^ * fe (r)a fe exp 



(A3) 



Any two point correlation function of products of two non-condensate fields in Eqs. (|A1[) and (|A2I) can thus be written 
in terms of the average occupation of different non-condensate single particle states \^k) & 3~±, e.g.: 



*^(r,r)#x(r ,0) 



(N-No) 



J2 n(^k(r)(N k )(N - N ,T) exp 



ie k T 



(A4) 



.(N-No) I , \{N-No) 

where we used that (a[a; )_ = (^a J k akJ Ski = (Nk)(N — N ,T)Ski ■ The function 



(N k )(N -N ,T) = Tr^ \ a\a k Q N _ No 



-PHa 



2±(N-N ) 



Qn-No 



(A5) 



describes the average many particle occupation of a non- 
condensate single particle state \^k), given that (N~Nq) 
particles are in the non-condensate, and given a temper- 
ature T of the external heat reservoir. For the explicit 
derivation of analytical expressions for the occupation 



numbers (Nk)(N — Nq,T), see appendix 151 

Anti-normally ordered products of two point correla- 
tion functions of two non-condensate fields in the interac- 
tion picture arising in Eqs. (|A1[) and (| A2|) can be obtained 
correspondingly, turning into 



*x(r, r)*i_(f,0) 



(N —Nq) 



= £ ^k(r)n(r ) [{N k )(N — Nq,T) + 1] exp 
fc/o 



ie k T 



(A6) 



where we have used that 



/ +\(N-N ) 

(akaj) = [(N k )(N-N ,T) + l]S kl 



(A7) 



Hence, we find for normally and anti-normally ordered two point correlation functions with respect to single particle 
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processes: 



g^(v,r',N-N ,T,T)=2 £ tJ(r^ t (r>[(f)*,(r>;(rX(0 [{N k )(N - N Q , T) + 1] 

~i(e k -e t - e m )r~ 



x (Ni)(N - N ,T){N m )(N - N ,T) exp 
+ 4 ^ |* fc (r)| 2 *r(r)vI/ i (f')|* m (f')| 2 (7V fc )(7V-7Vo,r) ; 

k,l,m^0 



x (JV I )(JV-JV 0j T)(JV m )(JV-JV ,T)ejq) 



(A8) 



g(-\r,v',N-N ,T,r)=2 ]T *J(?)* Jk (?')*f(?)* I (i ! ')*^(f)* m (?')(JV fc >(JV - JV 0j T): 

fc,Z,m^0 

x [(JV,)(JV - JV ,T) + 1] [(N m )(N - JV ,T) + 1] exp [ - {ek - % - £m) T 
+ 4 £ |* fc (r)| 2 *f(r)* ; (f')|vI/ m (r)| 2 (7V fc )(^-^o,r)x 

k,l,m=£0 



(A9) 



x [(JV,)(JV-JV ,r) + l] (N m )(N-N ,T) exp 

I 

Integration of (r, f , A^-iV , T, r) over the time inter- Next, we decompose the correlation functions for pair 

val t, multiplied by ^(r^j^r )exp[±i/ior], which arises events, Q<£) (r, r , N - N , T, r). Using Eq. dH EH, the 

from the backswitch of the condensate fields from the normally ordered correlation function for pair events is 

interaction picture, leads to the single particle loss and given by: 
feeding rates in Eq. ([57]) . 



g£)(r,r,N-N ,T,T) = r)*t(F, r)*_ L (r , 0)£ x (f , 0) 



(iV-JV ) 



/ - + \ (JV-iVo) / „ . „ \ (AT-JVo) 

2(^(r,T)*x(r,0)) (^(r, T)* ± (f, 0) ) 

2 n(^k(r)^n^i(r)(N k )(N - N ,T)(N l )(N - N Q ,T) exp 



+i (e k + ei) t 



The anti- normally ordered pair correlation function Q^} (f,r,N — No,T, r) can be decomposed similarly: 



(A10) 



Gt) (r, ?,N- N , T, t) = (* x (?, t)* x (?, t)^^, 0)*+_(f , 0) 



(N-N ) 



2(* ± (?,t)*5.(?',0)) (* x (?,r)*5_(?',0)) 

2 ]T ^(^(^(^v^ exp 



(All) 



which again, after multiplication with ^(r) (^(r')) exp[±2i/zor] and integration over r, turns into the pair feeding 
and loss rates in Eq. (|45|). 



Appendix B: Single particle non-condensate 
occupations 



(N k )(N - Nq,T) in Eq. (|A"5|) . for each particular non- 



The state of the non-condensate in Eq. (|22|) allows 
to determine the average number of particles, (N k ) = 
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condensate single particle mode \^f k ), given that Nq par- 
ticles populate the condensate mode, and consequently 
(TV — No) particles populate the non-condensate single 
particle modes. According to Eq. (|A5|) . we hence con- 



sider the expectation value of the number operator TVfc 
in a non-condensate state of (-/V — Nq) particles, which 
leads to 



(N-N ) 

(N k )(N-N ,T) = Z± 1 (N-N ) ^ ex P 

{N k } 



k^O 



(Bl) 



where Z±(N — No) is the partition function of (N — 
No) indistinguishable particles in the non-condensate in 
Eq. ([24| . In terms of the partial partition sum, Z± (TV — 



No) [25j|, which excludes the sum over one particular non- 
condensate single particle mode \^ k ), Eq. (|B1[) can be 
written as 



(N-N ) 



(N k )(N~No,T) = Z- 1 (N-N ) ^ N k exp [-0e k N k ] zf>{N - N - N k ) 



(B2) 



N k =0 



For small enough N k (it suffices to start at N k = 1 and to determine Z±\n — No — Nk) iteratively) , we can expand 



In 



Z { 1 ] (TV - N - 1) ~ln Zf (N - No) - (N - N ,T) , 



r(fe) 



(fe)/ 



which introduces the parameter 

a ( f>(N-N 0> T) 

From Eq. (|B3j) . we hence find the recursion relation 

z[ fc) (TV-TV - 1) 



din 



zf\N - No) 



Z { l\N -No) 



exp 



8(N - No) 



a W(N-N ,T) 



(B3) 



(B4) 



(B5) 



between the partial partition sums Z^ (N — No) of (TV — TVo) and Z± (TV — TVo — 1) of (TV — TVo — 1) non-condensate 
particles. Multiple iteration of Eq. ()B5|) hence leads to 

z[ fc) (TV-TV -TV fc 



z| fe) (TV- TV ) 



and Eq. (|B2[) turns into 



(TV fc )(TV-TV ,T) = 



Z^(N-N ) 
Z ± (N-N ) 



exp 



(N-N ) 



-N k a<£\N-No,T) 



(B6) 



N k =0 



J2 Nkcxp -((3e k +a^(N-N ,T)jN k 



(B7) 
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It remains to apply the same procedure to the partition function Z±(N — No). Using the decomposition in Eq. (|B7|) 
and applying Eq. (|B6[) . one finds that 



(B8) 



Setting Eq. (|B8|) into Eq. (|B1|) . the expectation value of particle number occupations of a particular non-condensate 
single particle state \^k), given that (N — No) particles are in the non-condensate, is given by 



(N-N ) 

Z±(N - N ) = Z [ l ] (N - N ) exp[-(p€ k + a^\N-No,T)^N k 

N k =0 



(N k )(N-N ,T) = 



1 



cxp 



Pe k + a { V(N-N Q ,T) 
I 



(B9) 



We now use that the parameter a±\N — Nq , T) is ap- 
proximately independent of the state k [231, 1 - e - tnc 
change of non-condensate single particle number occupa- 
tions during condensation is described by one single pa- 



J 



rameter, ~ a±(N — Nq,T), which is determined by 
the constraint of particle number conservation, as spelled 
out by the implicit equation 



J2(N k )(N-No,T) = J2 



fc^o 



exp[Pe k + a ± (N - N 0l T)} _ 1 



(N-N ) 



(BIO) 



r 



As evident from Eq. (|B4[) , and the fact that each subspace 
of (N — No) particles is in a thermal state, the parameter 
aj_(N — N ,T) can be interpreted as the ratio of the 
non-condensate chemical potential for a state of (N — 
N ) atoms, to the thermal energy Hence, from the 
definition in Eq. (|B4[) , we see that a±(N — N ,T) is, 
upon a constant, nothing more than the derivative of the 
Helmholtz free energy, a±(N-N ,T) = -/3~HnZ±(N - 



iVo) of the (N — Nq) particles in the non-condensate |25| . 
related to n ± (N — Nq,T) by 



a ± (N - N ,T) = -Phj_(N - N ,T) , 



(Bll) 



which introduces the non-condensate chemical potential 
ft±{N-N ,T). 
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